
clc; clear; close all;

%% 参数设置
N = 500000;  % 迭代步数
% 各模型采用的频率数组，可根据需要调整
omega_QIM = [0.80, 1.50, 1.00, 1.80];
omega_SIM = [0.50, 1.50, 1.00, 1.80];
omega_TGM = [0.10, 0.80, 0.40, 1.00];
omega_CGM = [0.10, 0.20, 0.15, 0.25];

% 颜色列表
color_list = {'r','g','b','k','m','c'};

%% Q-IM 模型绘图
subplot(2,2,1);
hold on; grid on; box on;
legend_str = cell(length(omega_QIM),1);
for idx = 1:length(omega_QIM)
    w = omega_QIM(idx);
    A = 1; % 幅值
    [v, i] = QIM(w, N, A);
    plot(v, i, color_list{idx}, 'LineWidth',1);
    legend_str{idx} = sprintf('\\omega=%.2f', w);
end
xlabel('v_1(n), V'); ylabel('i_1(n), A');
title('(a) Q-IM');
legend(legend_str, 'Location','best', 'NumColumns',2);

%% S-IM 模型绘图
subplot(2,2,2);
hold on; grid on; box on;
legend_str = cell(length(omega_SIM),1);
for idx = 1:length(omega_SIM)
    w = omega_SIM(idx);
    A = 0.1; % 幅值
    [v, i] = SIM(w, N, A);
    plot(v, i, color_list{idx}, 'LineWidth',1);
    legend_str{idx} = sprintf('\\omega=%.2f', w);
end
xlabel('v_2(n), V'); ylabel('i_2(n), A');
title('(b) S-IM');
legend(legend_str, 'Location','best', 'NumColumns',2);

%% T-GM 模型绘图
subplot(2,2,3);
hold on; grid on; box on;
legend_str = cell(length(omega_TGM),1);
for idx = 1:length(omega_TGM)
    w = omega_TGM(idx);
    A = 1;
    [v, i] = TGM(w, N, A);
    plot(v, i, color_list{idx}, 'LineWidth',1);
    legend_str{idx} = sprintf('\\omega=%.2f', w);
end
xlabel('v_3(n), V'); ylabel('i_3(n), A');
title('(c) T-GM');
legend(legend_str, 'Location','best', 'NumColumns',2);

%% C-GM 模型绘图
subplot(2,2,4);
hold on; grid on; box on;
legend_str = cell(length(omega_CGM),1);
for idx = 1:length(omega_CGM)
    w = omega_CGM(idx);
    A = 1;
    [v, i] = CGM(w, N, A);
    plot(v, i, color_list{idx}, 'LineWidth',1);
    legend_str{idx} = sprintf('\\omega=%.2f', w);
end
xlabel('v_4(n), V'); ylabel('i_4(n), A');
title('(d) C-GM');
legend(legend_str, 'Location','best', 'NumColumns',2);

sgtitle('四种离散忆阻器的钳制回滞曲线');

%% 封装函数：Q-IM 模型
function [v_array, i_array] = QIM(omega, N, A)
    % 初始状态
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = (phi^2 - 2) * v;
        v_array(n) = v;
        i_array(n) = i;
        % 更新 phi
        phi = phi + v;
    end
end

%% 封装函数：S-IM 模型
function [v_array, i_array] = SIM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = sin(pi * phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        phi = phi + v;
    end
end

%% 封装函数：T-GM 模型
function [v_array, i_array] = TGM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = sin(phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        % phi 更新：带有非线性项
        phi =  0.12 * phi * abs(phi) + 0.86 * phi - 0.02 * (phi^3) + 0.1 * v;
    end
end

%% 封装函数：C-GM 模型
function [v_array, i_array] = CGM(omega, N, A)
    phi = 0;
    v_array = zeros(1, N);
    i_array = zeros(1, N);
    for n = 1:N
        v = A * sin(omega * n);
        i = cos(phi) * v;
        v_array(n) = v;
        i_array(n) = i;
        % phi 更新：带 tanh 非线性项
        phi = 0.1 * tanh(phi^3) + 0.99 * phi + 0.1 * v;
    end
end
